######################################################
# South Korea Survey Experiment #
######################################################
# set working directory 
library(survival); library(psych); library(tidyverse);
library(effects);library(plyr)
library(stargazer);library(xtable); library(haven);
library(texreg);library(rstatix);library(MASS)

# load the data 
data <- read_sav("data.sav")
#----------------------------------------------------------#
# pre-treatment variables #
#----------------------------------------------------------#
# female
colnames(data)[5] <- "female"
data$female <- ifelse(data$female == 2, 1, 0)
# age
colnames(data)[6] <- "age"
data$age[data$age > 65] <- NA
data$age <- scale(data$age)
# education
colnames(data)[7] <- "education"
# college
data$college <- ifelse(data$education > 2, 1, 0)
data$education <- scale(data$education)
# income
colnames(data)[52] <- "income"
data$income <- scale(data$income)
#----------------------------------------------------------#
# Treatment # 
#----------------------------------------------------------#
# treatment: 1 control / 2 automation / 3 globalization
data$treatment <- ifelse(data$rnd == 1, 0, 1)
data$automation <- ifelse(data$rnd == 2, 1, 0)
data$globalization <- ifelse(data$rnd == 3, 1, 0)
# control
data$control <- NA
data$control <- ifelse(data$rnd == 1, 1, 0)
#----------------------------------------------------------#
# dependent variables
#----------------------------------------------------------#
data$E15[data$E15 == 9] <- NA
data$consumption <- 5 - data$E15
# scale 
data$con <- scale(data$consumption)
# recoding for binary
data$Rcon <- ifelse(data$E15 < 3, 1, 0)
write.csv(data, "data.csv")
# read new data 
data <- read.csv("data.csv")
###############################
# Figure 3a # 
###############################
# data for figure 3a
fig3a <- emmeans_test(data = data, formula = Rcon ~ treatment, p.adjust.method = "fdr")
fig3a <- get_emmeans(fig3a)
fig3a$treatment <- ifelse(fig3a$treatment == 1, 'Treatment', 'Control')
# plot for figure 3a
ggplot(fig3a, aes(x = treatment, y = emmean, fill = treatment)) +
  geom_bar(stat = "identity") +
  ylab("Proportion of respondents \n who agree or lean to social consumption") +
  xlab("") +
  ylim(0,1) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0.2, position = position_dodge(18)) +
  scale_fill_manual(values=c("#999999", "#32B141")) +
  guides(fill = "none") +
  theme_light() +
  theme(axis.title.y = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16)) +
  geom_text(
    aes(x = treatment, y = emmean, label = round(emmean, 2)),
    position = position_stack(vjust = 0.5),
    size = 4)
ggsave("figure3a.pdf")
###############################
# Figure 3b # 
###############################
# data for figure 3b
fig3b <- emmeans_test(data = data, formula = Rcon ~ rnd, p.adjust.method = "fdr")
fig3b <- get_emmeans(fig3b)
fig3b$treatment <- NA
fig3b$treatment[fig3b$rnd == 1] <- "Control"
fig3b$treatment[fig3b$rnd == 2] <- "Automation"
fig3b$treatment[fig3b$rnd == 3] <- "Globalization"
# plot for figure 3b
ggplot(fig3b, aes(x = treatment, y = emmean, fill = treatment)) +
  geom_bar(stat = "identity") +
  ylab("Proportion of respondents \n who agree or lean to social consumption") +
  xlab("") +
  ylim(0,1) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                width = 0.2, position = position_dodge(18)) +
  scale_fill_manual(values=c("#72A050", "grey", "#32B141")) +
  guides(fill = "none") +
  theme_light() +
  theme(axis.title.y = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16)) +
  geom_text(
    aes(x = treatment, y = emmean, label = round(emmean, 2)),
    position = position_stack(vjust = 0.5),
    size = 4)
ggsave("figure3b.pdf")
###############################
## Table 2  ##
###############################
m1 <- lm(con ~ treatment, data = data)
m2 <- lm(con ~ treatment + education + income +  
           + female + age, data = data)
m3 <- lm(con ~ globalization + automation, data = data)
m4 <- lm(con ~ globalization + automation + education + income +  
           + female + age, data = data)
texreg(list(m1, m2, m3, m4))
###############################
## Appendix Table A7&8 ##
###############################
# Table A.6 summary statistics
survey.demo <- subset(data, select = c(age, female, education, income, consumption))
survey.demo <- describe(survey.demo)
survey.demo <- subset(survey.demo, select = c(n, mean, median, sd, min, max))
xtable(survey.demo)
###############################
## Appendix Table A9 ##
###############################
# balance test 
t.test(data$age[data$treatment == 1], data$age[data$treatment == 0])
t.test(data$female[data$treatment == 1], data$female[data$treatment == 0])
t.test(data$income[data$treatment == 1], data$income[data$treatment == 0])
t.test(data$education[data$treatment == 1], data$education[data$treatment == 0])
###############################
## Appendix Table A10 ##
###############################
dt <- read_sav("data.sav")
colnames(dt)[5] <- "female"
dt$female <- ifelse(dt$female == 2, 1, 0)
colnames(dt)[6] <- "age"
dt$age[dt$age > 65] <- NA
colnames(dt)[7] <- "education"
colnames(dt)[52] <- "income"
dt$treatment <- ifelse(dt$rnd == 1, 0, 1)
dt$automation <- ifelse(dt$rnd == 2, 1, 0)
dt$globalization <- ifelse(dt$rnd == 3, 1, 0)
dt$E15[dt$E15 == 9] <- NA
dt$consumption <- 5 - dt$E15

a1 <- lm(consumption ~ treatment, data = dt)
a2 <- lm(consumption ~ treatment + education + income +  
           + female + age, data = dt)
a3 <- lm(consumption ~ globalization + automation, data = dt)
a4 <- lm(consumption ~ globalization + automation + education + income +  
           + female + age, data = dt)
texreg(list(a1, a2, a3, a4), style = "AJPS")
###############################
## Appendix Table A11 ##
###############################
# ordered probit 
data$Fcon <- as.factor(data$consumption)
p1 <- polr(Fcon ~ treatment, data = data)
p2 <- polr(Fcon ~ treatment + education + income + female + age, 
     data = data)
p3 <- polr(Fcon ~ globalization + automation, data = data)
p4 <- polr(Fcon ~ globalization + automation + education + income +  
       + female + age, data = data)
texreg(list(p1, p2, p3, p4))
###############################
## Appendix Table A12 ##
###############################
# placebo test 
data$E12_4[data$E12_4 == 9] <- NA
data$NE12_4 <- NA
data$NE12_4 = (data$E12_4 - 1)/3
data$E12_5[data$E12_5 == 9] <- NA
data$NE12_5 <- NA
data$NE12_5 = (data$E12_5 - 1)/4

p1 <- lm(NE12_4 ~ treatment, data = data)
p2 <- lm(NE12_4 ~ treatment + education + income +  
           + female + age, data = data)
p3 <- lm(NE12_5 ~ treatment, data = data)
p4 <- lm(NE12_5 ~ treatment + education + income +  
           + female + age, data = data)
texreg(list(p1, p2, p3, p4), style = "AJPS")
###############################
## Appendix Table A13 ##
###############################
c1 <- lm(con ~ globalization + automation, data = data[data$college == 1,])
c2 <- lm(con ~ globalization + automation + income +  
           + female + age, data = data[data$college == 1,])
c3 <- lm(con ~ globalization + automation, data = data[data$college == 0,])
c4 <- lm(con ~ globalization + automation + income +  
           + female + age, data = data[data$college == 0,])
texreg(list(c1, c2, c3, c4), style = "AJPS")

con_edu <- data.frame(
  Education = factor(rep(c("College or Above", "Below College"), each=2)),
  estimate = c(0.18, 0.08, 0.34, 0.07),
  sd = c(0.11, 0.11, 0.14, 0.14),
  Type = rep(c("Globalization", "Automation"), 2)
)

plot1 <- ggplot(con_edu, aes(y = estimate, x = Education, group = Type)) +
  geom_point(position =  position_dodge(width = 0.5)) +
  geom_linerange(aes(ymin = estimate-1.96*sd, ymax = estimate+1.96*sd),
                 position = position_dodge(width = 0.5), lwd = 1/2) +
  geom_linerange(aes(ymin = estimate-1.645*sd, ymax = estimate+1.645*sd),
                 position = position_dodge(width = 0.5), lwd = 1) +
  labs(x = "", y = "Coefficient estimate", color = "") +
  geom_hline(yintercept = 0,linetype = "dashed") +
  theme_bw(base_size = 14)+
  ylim(-0.6, 0.7) + 
  geom_text(aes(label = Type), position = position_dodge(width = 0.5), vjust = -1, hjust = -0.1)
ggsave("con_edu.pdf", width = 9,  height = 7, dpi = 1200)

###############################
## Appendix Table A14 ##
###############################
data$highinc <- ifelse(data$income > -0.1821749, 1, 0)
c1 <- lm(con ~ globalization + automation, data = data[data$highinc == 1,])
c2 <- lm(con ~ globalization + automation + education +  
           + female + age, data = data[data$highinc == 1,])
c3 <- lm(con ~ globalization + automation, data = data[data$highinc == 0,])
c4 <- lm(con ~ globalization + automation + education +  
           + female + age, data = data[data$highinc == 0,])
texreg(list(c1, c2, c3, c4), style = "AJPS")

con_inc <- data.frame(
  Income = factor(rep(c("Above median", "Below median"), each=2)),
  estimate = c(-0.14, -0.21, 0.54, 0.30),
  sd = c(0.12, 0.12, 0.12, 0.12),
  Type = rep(c("Globalization", "Automation"), 2)
)

plot2 <- ggplot(con_inc, aes(y = estimate, x = Income, group = Type)) +
  geom_point(position =  position_dodge(width = 0.5)) +
  geom_linerange(aes(ymin = estimate-1.96*sd, ymax = estimate+1.96*sd),
                 position = position_dodge(width = 0.5), lwd = 1/2) +
  geom_linerange(aes(ymin = estimate-1.645*sd, ymax = estimate+1.645*sd),
                 position = position_dodge(width = 0.5), lwd = 1) +
  labs(x = "", y = "Coefficient estimate", color = "") +
  geom_hline(yintercept = 0,linetype = "dashed") +
  theme_bw(base_size = 14)+
  ylim(-0.8, 0.8) +
  geom_text(aes(label = Type), position = position_dodge(width = 0.5), vjust = -1, hjust = -0.1)
ggsave("con_inc.pdf", width = 9,  height = 7, dpi = 1200)

###############################
## Appendix Table A15 ##
###############################
data$unemp <- ifelse(data$E8 == 1, 1, 0)
c1 <- lm(con ~ globalization + automation, data = data[data$unemp == 1,])
c2 <- lm(con ~ globalization + automation + income + education + 
           + female + age, data = data[data$unemp == 1,])
c3 <- lm(con ~ globalization + automation, data = data[data$unemp == 0,])
c4 <- lm(con ~ globalization + automation + income + education + 
           + female + age, data = data[data$unemp == 0,])
texreg(list(c1, c2, c3, c4), style = "AJPS")

con_uemp <- data.frame(
  Unemployment = factor(rep(c("Unemployment experience (yes)", "Unemployment experience (no)"), each=2)),
  estimate = c(0.32, 0.22, 0.17, -0.05),
  sd = c(0.13, 0.13, 0.12, 0.12),
  Type = rep(c("Globalization", "Automation"), 2)
)

plot3 <- ggplot(con_uemp, aes(y = estimate, x = Unemployment, group = Type)) +
  geom_point(position =  position_dodge(width = 0.5)) +
  geom_linerange(aes(ymin = estimate-1.96*sd, ymax = estimate+1.96*sd),
                 position = position_dodge(width = 0.5), lwd = 1/2) +
  geom_linerange(aes(ymin = estimate-1.645*sd, ymax = estimate+1.645*sd),
                 position = position_dodge(width = 0.5), lwd = 1) +
  labs(x = "", y = "Coefficient estimate", color = "") +
  geom_hline(yintercept = 0,linetype = "dashed") +
  theme_bw(base_size = 14)+
  ylim(-0.6, 0.6) +
  geom_text(aes(label = Type), position = position_dodge(width = 0.5), vjust = -1, hjust = -0.1)
ggsave("con_uemp.pdf", width = 9,  height = 7, dpi = 1200)

